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Abstract 

The motion of a Cqq molecule over a graphene sheet at finite temperature is investigated 
both theoretically and computationally. We show that a graphene sheet generates a van der 
Waals laterally periodic potential, which directly influences the motion of external objects in its 
proximity. The translational motion of a Cqq molecule near a graphene sheet is found to be diffusive 
in the lateral directions. While, in the perpendicular direction, the motion may be described as 
diffusion in an effective harmonic potential which is determined from the distribution function of 
the position of the Cqo molecule. We also examine the rotational diffusion of Cqq and show that 
its motion over the graphene sheet is not a rolling motion. 
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1 Introduction 



Various properties of graphene as a new two dimensional material have been studied both experimen- 
tally and theoretically [1, 2]. A recent experimental research based on TEM visualization studies the 
images and dynamics of light atoms deposited on a single-layer graphene sheet [3] . On the other hand 
diffusion process and crystallization of atoms and light molecules on various surfaces have also been 
subject of research for many years due to both their theoretical importance as well as their techno- 
logical applications. The theoretical studies of the motion of molecular scale objects on the various 
surface is also applicable to the motion of nanoscale object's over nanoelectromechanical surfaces, i.e., 
graphene. These studies play an important role in designing graphene based nanosensors. Diffusive 
motion of inclusions (e.g. macromolecules) over a rough membrane is another related subject, which 
has received a lot of attention in recent years [4, 5, 6]. Solidification of Cqq molecules over various 
substrates particulary graphene is another subject of interest. The van der Waals epitaxy of a solid 
6*60 over graphene sheet has been done in recent experiments [7]. Also various dynamical properties 
of the spinning motion of 6*60 on the Au(l 1 1) surface have been studied by Teobaldi et al. by means 
of molecular dynamics (MD) simulations [8]. 

The choice of van der Waals parameters of Cgo-graphene interaction is an important issue for 
performing any molecular dynamics simulations. The C6o~carbon nanotubes interaction can be ob- 
tained both experimentally and theoretically [9]. Several van der Waals parameters for physical 
adsorption of C^o on graphite and other substrates were formulated using a continuum rigid body 
model for Ceo and a continuum dielectric media for graphite by Girad et al [10]. The charge transfer 
from graphene to the 6*60 molecule is an open question which may be tackled using modern density 
functional theories. For alkali metals on the graphite there have been some calculations in order to 
estimate transferred charges [11]. 

Two main approaches for investigating the diffusive motion of an external inclusion on the mem- 
branes surface can be introduced. The first one is based on the continuum description for the mem- 
brane structure and introduces an effective Hamiltonian which can be used to study curvature-coupled 
diffusion [5, 4, 6]. The elastic aspects of the membrane thus play a key role in this approach. The 
second approach is an atomistic one which considers the details of the membrane and inclusion struc- 
tures as well as atomic interactions on different levels of coarse graining [12, 13]. The main goal in 
these studies is to characterize the inclusion's motion. 

The diffusive motion of particles may be modeled generally via a Langevin equation. Lacasta et 
al have modeled general two dimensional solid surfaces using a two dimensional potential both with 
a deterministic periodic potential and a random one [14]. Different values for dissipation parameters 
generate different trajectories for a point-like particle over such potentials, which thus lead to different 
dynamical behaviors ranging from subdiffusive to superdiffusive motion [14]. 

In a previous study we showed that the alkali and transition metals distribute over graphene and 
bind to its surface via a Lennard- Jones potential and construct atomic nanoclusters [16]. For metallic 
nanoclusters, no diffusive motion was found in low temperature. It is also found that potassium 
atoms, in low temperature, construct a particular phase on the graphene sheet as well as on the 
graphite [11]. To our knowledge there have been only a few number of studies on the interaction 
between graphene and inclusions that address the pattern of inclusion motion [4, 5, 6]. 

In this paper we study the dynamics of a single molecule, specifically chosen as the Cqq molecule, 
on a graphene surface. We show that a graphene sheet creates a periodic van der Waals potential 
in its surrounding space and that there are some simple criteria to determine when a molecule may 
diffuse through the potential wells generated by the graphene (as assisted by thermal fluctuations) 
or be likely to be trapped in the potential wells near the surface. We introduce a laterally-averaged 
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effective potential for the graphene sheet from the distribution function of Cqq near the surface, and 
show that this potential may be approximated best by a harmonic potential in the direction normal 
to the sheet. We also introduce an effective friction coefficient for the diffusive motion of Cgo over 
graphene. Finally, we show that the motion of Cqq over graphene is not a rolling motion and also the 
variation of the perpendicular component of angular velocity of Cqo is greater than its component 
parallel to the graphene sheet. 

2 Methods 

We employ classical Molecular Dynamics (MD) algorithm to simulate the Cgo-graphene system. The 
graphene is modeled as a square-like sheet of area 1883 nm 2 constructed of 72000 carbon atoms. The 
temperature was kept constant (300K) in our simulations by employing a Nose-Hoover thermostat. 
For the covalent bounds between the nearest neighbor atoms of the graphene sheet (and for chemically 
bonded atoms in Cqo molecule) we have used Brenner's potential [17]. For the interaction between 
each atom of graphene with each of Cqqs atoms, we have used the Lennard-Jones potential (LJ), 

U LJ = 4e{(a/r) 12 -(a/rf} 

with a typical values a = 3.4 A and e = 2.4 meV [15]. A graphene sheet was initially positioned, on 
average, at z = plane, and the centers of the Cqq molecule were placed above it at z = 7 A. The 
temperature-dependent molecular dynamics simulations run for up to two nanoseconds. 




Figure 1: (Color online) Two common triangular sublattices of graphene sheet in the x — y plane. 



3 Two dimensional potential of garphene sheet 

First we study the periodic two dimensional LJ potential that a flat graphene sheet creates near its 
surface. Generally, graphene is not a flat sheet in finite temperature and exhibits small roughness; it 
may thus be locally approximated with a flat sheet [18]. The aforementioned quantity helps one to 
gain some insight into how graphene influences molecules such as Cqq moving at its proximity. 

A flat graphene sheet comprises two triangular Bravais sublattices (see Fig. 1), and the LJ potential 
due to this sheet [9] can be written as the sum of these two sublattice potentials 

E T (x, y, z) = VA-lattUx + VB-lattice, (1) 
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Figure 2: (Color online) Periodic two dimensional potential, V(x,y), created by the graphene sheet. 
Here x and y refer to spatial position at the height zq=3.5A above the graphene sheet. This figure 
shows the top view of the potential. 



where (x, y, z) is the position for center of mass of the molecule, having N atoms, above the sheet. 
We may write the appropriate expression for VA-iattice as 

y _ y^y^y^ ^ ' (2) 

A-latUce ((x + xi _ n _ m/2)2 + {y + y[ _ ^ m/2 )2 + {z + Z{ )2)6/fc 

and for Vs-iattice as 

y _ y^y^y^ ^ > 

B ~ latUCe knhtl « x +Xi + l/2-n- m/2)2 + {y + m - V3m/2 - V3/6) 2 + (z + z/) 2 ) 6 / fe 

Here m and n are integer numbers which count lattice points for each sublattice. The coordinate 
of the I th atom (xi,yi,z{) is measured from the molecule's center of mass. The sum over k also, is 
responsible for varying of fraction's power between 12 or 6, and switching its sign. Because of the 
short range behavior of the LJ potential, we observed that using the cut-offs of \m\ > 10 or \n\ > 10 
leads to quite accurate results with negligible cut-off errors in the total potential/force. For simplicity 
in Eqs. (2) and (3) all lengths were re-scaled by a\ = \/3ao where ao = 1.42 A and energies are in 
units of e, with a\ being the length of primitive vector of sublattices (see Fig. 1). 



3.1 Potential energy between a single carbon atom and graphene sheet 

At a fixed height above the graphene sheet, z = zq, the total potential in Eq. (1) reduces to a two- 
dimensional potential V(x,y), which is periodic in x and y. Assume that we have fixed the height 
of a point like particle, such as a carbon atom, at a given value. Figure 2 shows this periodic two- 
dimensional potential for zq=3.5 A. This potential provides some insights into the diffusive motion 
or trapping of a point-like particle close to the graphene sheet [3]. The required potential can be 
obtained from Eqs. (2) and (3) with N = 1 and putting xi, y\ and z\ equal to zero. Figure 3 shows 
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Figure 3: (Color online) Periodic potential energy between flat graphene sheet and a single carbon 
atom at x=0 as a function of y for several height values. 



the variation of potential energy in y direction at x = for five different heights. To show that 
the variation of V(x, y) in the z direction appears to be similar to the LJ potential, with different 
functionality, we averaged the above potential in the Wigner-Seitz primitive cells on both sublattices 
and plot the result in the bottom panel of Fig. 4. As can be seen from the figure the potential 
minimum between graphene sheet and the single carbon atom is deeper than the simple case of two 
interacting carbon atoms via the LJ interaction. Furthermore, the minimum distance (z m i n ~3.4 A) 
is not z = 2 1 / 6 <7 ~ 3.8 A as it is in the usual LJ (solid blue curve in Fig. 4). The minimum point of 
V z gives the equilibrium distance of a carbon atom over graphene flat sheet. 

Note that the particle will equilibrate with the graphene sheet and obtains a mean kinetic energy 
of the order of KbT in its 2D lateral motion. If KbT ~ 25.7 meV is smaller than the potential barriers 
heights, the particle may be trapped in one of the potential wells. Since V(x, y) is a two-dimensional 
potential, the particle could take various paths from one well to neighboring wells. But for simplicity 
we may roughly take the difference between maximum and minimum values of V(x, y) appearing in 
Wigner-Seitz primitive cells as a measure of the barrier strength compared to the thermal energy 
KbT. The top panel in Fig. 4 shows the difference between maximum and minimum AV in units of 
KbT for various values of the particle height. 

3.2 Potential energy between 6*60 molecule and graphene sheet 

In the case of molecules such as 6*60 as considered here, the total two-dimensional potential is periodic 
as well, see Fig. 5. The interaction between this molecule and the graphene sheet has several functional 
forms depending on the orientation of the molecule over the sheet even at a fixed height of the 
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Figure 4: (Color online) Top: Difference between maximum and minimum of two-dimensional po- 
tential (averaged over Wigner-Seitz primitive cell of both A-lattice and B-lattice) created by a flat 
graphene sheet (see Fig. 2) as a function of normal distance z as experienced by a single carbon atom. 
Bottom: Total potential energy between a flat graphene sheet and a single carbon atom has been 
depicted. The solid blue curve is a simple LJ potential between two carbon atoms. In both panels, to 
eliminate the dependence on the x and y variables, we averaged over the first Brillouin zone of both 
A-lattice and B-lattice. The data have small (10 -5 ) error bars which are not shown. 



molecule's center of mass above the sheet. Two particular orientations of Cqq are more interesting 
than the others. These two refer to the cases when a pentagon or a hexagon of Cqo's atoms faces the 
flat graphene surface. Figure 5 is related to a Cqo molecule when one of the pentagons is near to the 
surface and the plane of the pentagon is parallel to the graphene sheet. Since the radius of the Cqq 
molecule is Rc 60 =3.54 A, it can never get closer to the surface than this distance. Figure 6 shows 
the variation of potential energy along normal direction for five different height of the center of 
mass of 6*60 molecule. Furthermore the variation of Et along z direction averaged over Wigner-Seitz 
primitive cells of A-lattice and B-lattice is shown in bottom panel of Fig. 7. This gives the minimum 
height value for the center of mass of Cqq molecule as z m i n ~ 6.5 A which is obviously larger than 
the equilibrium distance obtained for a single carbon above the flat sheet. The binding energy of 
Cqo molecule and monolayer graphene can be estimated around 800 mev which is the same to the 
experimental value for binding energy of Cqq molecule and graphite bulk [9] . 

Similar to the the single carbon case we show the difference between maximum and minimum 
appearing in Wigner-Seitz primitive cells with respect to the thermal energy in top panel of Fig. 7. 
The several orientations have not much effects in the curves in Fig. 7. Furthermore, the height of 
saddle points in the potential profile in Fig. 5 with respect to the thermal energy is around 0.2, which 
indicates that thermal energy is the dominant factor. Therefore, we expect that there should not be 
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Figure 5: (Color online) Periodic two dimensional potential created by the graphene sheet. In this 
case x and y are the center of mass coordinates of Cqo measured from an origin on the sheet at the 
height z = 6.5A. This is obtained by summing the potential experienced by all individual atoms 
comprising the Cqo molecule. The panel in the left shows the top view of the potential, and the right 
panel shows the side view. The middle panel shows the potential energy between Cqo molecule and 
flat graphene sheet at x=0 as a function of y for several height values. 



any trapping in the motion of Cqo over graphene even at low temperatures and various orientations. 
In the next section we show that the motion coincides with a diffusive Brownian motion. 

This foregoing observation about the size of Cqo molecule and the ratio between thermal energy 
and the depth of the potential wells may be extended to other macromolecules and nanostructures 
such as bucky-balls and carbon nanotubes. Almost all of these carbon allotropes have longitudinal 
and lateral dimensions larger than those of Cqo- In a future work we will investigate dynamics of 
other carbon allotropes over this mono layer sheet at finite temperatures. 

4 Trans lational diffusion 

In Fig. 8, the mean square displacement (R 2 ) of an ensemble of 30 Cqo molecules moving near a 
graphene sheet at room temperature is plotted as a function of time for the projected two dimensional 
motion onto the x-y plane, where R 2 = x 2 + y 2 . The inset of the figure shows x — y trajectory for a 
single Cqo molecule. The total simulation time is 2.5 ns. Diffusion coefficient for the Cqo molecule is 
obtained from this graph as D = 7.0 x 10 m 2 s _1 . Using Einstein's relation, we estimate an effective 
friction coefficient of £ = = 1.4 x 10 9 m~ 2 s x KbT for Cqo moving near a graphene sheet. At very 
short times the observed motion is not diffusive because we put all Cqo at the center of the membrane 
with initial random velocities extracted from a Maxwell-Boltzman distribution, and the graphene is 
not in thermal equilibrium with those molecules. Therefore at very short time of about 5 ps, the 
Ceo molecules attempt to find the minimum energy trajectories, yet their total displacement is not 
significant. 
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Figure 6: (Color online) Periodic potential energy between 6*60 molecule and flat grapheen sheet at 
x=0 as a function of y for several height values. 



5 Effective potential for vibrational motion in z-direction 

In our MD simulations the Cqo molecule does not unbind from the graphene sheet. We may introduce 
an effective potential for the Ceo motion in the z-direction by calculating the distribution function, 
p(z), of the height of a Cqq above the sheet. By virtue of the Boltzmann factor, we may define the 
effective potential as ^ T = — log[p(z)] which embeds in itself both entropic and energetic factors 
related to the equilibrated motion of Cgo near the graphene sheet. Results are represented in Fig. 9. 
It turns that the effective potential has a harmonic-type shape and ^ T = ^k(z — z m ) 2 with the 
effective parameters k and zq obtained from a best fit as k = 1.255-fTgT (N/A) and z m = 5.908 A. 
One thus expects that the Cqo molecules effectively exhibit bounded vibrational motion in the normal 
direction to the graphene sheet. In lateral directions, as we discussed in the previous sections, the 
motion is diffusive. Here the mean value for (z) =5.975 A shows the equilibrium distance of center 
of mass of Cqq molecule which already has been thermalized with the graphene sheet. This value of 
height depends on the temperature of the system. Note that this equilibrium distance is calculated 
for the rough graphene unlike the value reported in section 3. The distance from graphene sheet is 
thus found to be bigger than the distance of Cqq from the first layer of gold bulk [8], i.e. 5.4 A. This 
is because the graphene sheet is a single atomic layer instead of a bulk material. 

6 Rotational diffusion 

Rotation matrix transforms the coordinates of a vector in body-fixed frame to the coordinates in the 
lab-fixed frame. We define the body-fixed frame as two perpendicular vectors, which are chosen as two 
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Figure 7: (Color online) Top: Difference between maximum and minimum of two-dimensional poten- 
tial created by a flat graphene sheet (see Fig. 5) as a function of normal distance z as experienced 
by Cqo molecule . Bottom: Total potential energy between a flat graphene sheet and a Cqq molecule. 
In both panels, to eliminate the dependency of the x and y variables, we averaged over Wigner-Seitz 
primitive cells of both A-lattice and B-lattice with the data have infinitesimal (10 -5 ) error bars which 
were not shown. 



lines connecting two pairs of opposite points on the C§q cage and the cross vector of those two vectors. 
In order to be sure about the orthonormality of those vectors the Gram-Schmidt orthonormalization 
processes was applied in each time step of simulation. 

The quaternion representation of rotation matrix is a useful method for numerical simulations 
since they are more numerically stable. [19]. The angular velocities of Cqq in the body-fixed frame 
can be written in the term of quaternions and their time derivatives. Applying the inverse of the 
rotation matrix on the angular velocity of Cqo in the body-fixed frame, yields the angular velocities 
of 6*60 in lab- fixed frame, ~ut. 

To show that motion of 6*60 on the graphene is not a rolling motion, we have calculated the 
cross-correlation of the unit vector of angular velocity, 0, and the unit vector of velocity, v, i.e, 
{Q ■ v) = 3.0 x 1(T 2 ± 0.01 and {\u ■ v\) = 0.47 ± 0.01. Note that the first average is almost zero 
but it can not confirm that these two vectors are uncorrelated because this average becomes zero for 
two perpendicular vectors as well. On the other hand, the non-zero second average shows that they 
are not perpendicular. Therefore, these two averages together make sure that these two vectors are 
uncorrelated. The independence of the direction of velocity from the angular velocity ensures that 
the motion of Ceo over graphene sheet is not a rolling motion. 

Defining a fixed vector in Cqo, ji, helps us to investigate the diffusive nature of C§o motion . It 
can be understood from the time auto correlation of \i where the correlation time is 37.5 ps (see Fig. 
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Figure 8: (Color online) Mean square displacement, (R 2 ), versus time for the lateral motion of an 
ensemble of 30 Cqo molecules moving near a graphene sheet at room temperature. Inset show a 
typical x — y trajectory of a single Cqo on the graphene sheet. 



10). It is obvious that the correlation length of [i is about 75000 steps which means that after this 
time the orientation of Cqo becomes completely different. Angular velocity in different directions has 
different behaviors. Figure 11 shows the increment of lo z and 0J X y Obviously the fluctuations of lo z 
is bigger than uj xy . This is due to the structure of Cqo which is not a continues ball. It is a discrete 
spherical object. As we mentioned in the previous section the pairwise interaction of Cqo with the 
graphene atoms depends on the several orientations of Cqo- For example when a hexagon of Cqo is 
parallel to another hexagonal in the graphene plane, Cqo would be more stable than the other possible 
orientations. These restrictions do not affect the rotation around the z axis, so the rotation around 
the z axis is easier than around x — y direction. 

7 Conclusion 

We studied the motion of Cqo molecule over a graphene sheet. Both flat approximation for monolayer 
graphene sheet and monolayer graphene at a finite temperature have been studied using atomistic 
simulations. The depth of the potential wells generated by a graphene sheet in its proximity decreases 
as the height of an external object increases above the sheet. The binding energy of Cqo over a 
graphene sheet was found as 800 meV close to the experimental value for binding energy of Cqo 
and graphite [9]. The motion of Cqo in the perpendicular direction was found to be a vibrational 
motion similar to a simple harmonic oscillator. While the motion in lateral directions is found to be 
a diffusive non-rolling motion. 
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Figure 9: (Color online) Effective potential of the interaction between 6*60 molecule and a graphene 
sheet. 
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